Appendix 3

1 Overview

This document contains the code tutorial from “A comparison of statistical models used to characterize species-habitat associations with movement data”. This code tutorial demonstrates the analysis of ringed seal movement in eastern Hudson Bay, Canada.

2 Set-up

2.1 Load packages

library(tidyverse)
library(ggplot2)
library(lubridate)
library(rnaturalearth)
library(rgdal)
library(terra)
library(amt)
library(momentuHMM)
library(viridis)
library(sf)
library(here)

2.2 Fish data

Next we load in the fish data. This is a subset of the data from Florko et al. (2021).

# load fish data
fish <- read.csv(here("data/preydiv.csv"))

Visualize the fish data.

#prepare world data for mapping
## List of azimuthal equal area - HudsonBay 
natearth <- ne_countries(scale = "medium",returnclass = "sp")
natearth <- natearth[which(natearth$continent!="Antarctica"),]
nat_trans <- spTransform(natearth,CRS("+proj=laea +lat_0=60 +lon_0=-85 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"))

## Warning in spTransform(xSP, CRSobj, ...): NULL source CRS comment, falling back
## to PROJ string

## Warning in wkt(obj): CRS object has no comment

# plot fish map
fishmap <- ggplot() + 
  geom_tile(data = fish, aes(x = lon,y = lat,fill = preydiv)) +
  scale_fill_viridis(option = "mako", limits = c(0.5, 0.75), name = "Prey\ndiversity") +
  geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
  coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F) 
fishmap

2.3 Seal data

Next we load in the seal movement data. This is a subset of the data from Florko et al. (2023).

# load seal data
seal <- read.csv(here("data/seal_track_m.csv")) 
head(seal)

##        lon       lat             date       id
## 1 357940.9 -368949.6 2012-10-29 14:00 116484_1
## 2 371594.3 -353092.9 2012-10-30 14:00 116484_1
## 3 388429.2 -361478.4 2012-10-31 14:00 116484_1
## 4 399326.2 -357288.1 2012-11-01 14:00 116484_1
## 5 411788.2 -349500.7 2012-11-02 14:00 116484_1
## 6 407708.1 -372564.6 2012-11-03 14:00 116484_1

# ensure the data is in the correct format
seal <- seal %>%
  mutate(id = as.character(id),
         date = as.Date(date)) 

Visualize the seal data on top of the fish data.

# plot p
sealfishmap <- fishmap + 
  geom_point(data=seal, aes(x=lon, y=lat), alpha = 0.6, color = "#FCEEAE") + 
  geom_path(data=seal, aes(x=lon, y=lat), color = "#FCEEAE") 

sealfishmap

3 Fit models

We will fit four models: 1) a resource selection function (RSF), 2) a step selection function (SSF), 3) a SSF that allows the habitat covariate to affect movement, and 4) a hidden Markov model (HMM). All four of these models will include prey diversity as a covariate.

3.1 RSF

We will fit the RSF first, as it is the most simple of the four models. In preparation, we will generate an availability sample, create a map to visualize the spatial distribution of the used and available locations, and extract covariates for the used and available locations. Next, we will fit the model and view the model summary.

All of the step selection functions (both the RSF and two SSFs) will be fit using the ‘amt’ package (Signer et al. 2019).

# prep data and generate availability sample
set.seed(2023)
data_rsf <- seal %>%
  make_track(lon, lat, date) %>% # convert data to track format
  random_points() # generate availability sample; default is ten times as many available points as observed points

# plot used vs available locations on-top of prey diversity
sealfishmap +
  geom_point(data=data_rsf, aes(x=x_, y=y_, color = case_), alpha = 0.6) +
  scale_color_manual(values = c("black", "#FCEEAE"), 
                     label = c("Available", "Observed"), name = "Data type") 

# rasterize and extract prey diversity covariate
fish_raster <- terra::rast(fish, crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

data_rsf <- data_rsf %>%
  extract_covariates(fish_raster)

# fit rsf (a binomial logistic regression)
rsf1 <- data_rsf %>%
  amt::fit_rsf(case_ ~ preydiv, model = TRUE)

# see model summary
summary(rsf1)

## 
## Call:
## stats::glm(formula = formula, family = stats::binomial(link = "logit"), 
##     data = data, model = TRUE)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.5753  -0.4729  -0.4294  -0.3754   2.3363  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -6.118      1.400  -4.369 1.25e-05 ***
## preydiv        6.353      2.312   2.748    0.006 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 938.28  on 1539  degrees of freedom
## Residual deviance: 930.60  on 1538  degrees of freedom
## AIC: 934.6
## 
## Number of Fisher Scoring iterations: 5

3.2 SSF

Next we will fit our two SSFs. The workflow is similar to that of the RSF.

# prep data and generate availability sample
set.seed(2023)
data_ssf <- seal %>%
  make_track(lon, lat, date) %>% #  convert data to track format
  steps() %>% # convert track data to step format (i.e., with a start and an end)
  random_steps() %>% # generate availability sample 
  arrange(case_)

# plot used vs available locations on-top of prey diversity
sealfishmap +
  geom_point(data=data_ssf, aes(x=x2_, y=y2_, color = case_), alpha = 0.6) +
  scale_color_manual(values = c("black", "#FCEEAE"), 
                     label = c("Available", "Observed"), name = "Data type") 

# extract prey diversity covariate
data_ssf <- data_ssf %>%
  extract_covariates(fish_raster, where = "end") #sample at end of step

# transform movement covariates
data_ssf <- data_ssf %>%
  mutate(cos_ta = cos(ta_),
         log_sl = log(sl_))

# fit ssfs
## ssf1: ssf without covariate affecting movement kernel
ssf1 <- data_ssf %>%
  fit_clogit(case_ ~ preydiv + log_sl + cos_ta + strata(step_id_), model = TRUE)

## ssf2: ssf with covariate affecting movement kernel
ssf2 <- data_ssf %>%
  fit_clogit(case_ ~ preydiv*log_sl + cos_ta + strata(step_id_), model = TRUE)

# see model summaries
summary(ssf1)

## Call:
## coxph(formula = Surv(rep(1, 1518L), case_) ~ preydiv + log_sl + 
##     cos_ta + strata(step_id_), data = data, model = TRUE, method = "exact")
## 
##   n= 1516, number of events= 138 
##    (2 observations deleted due to missingness)
## 
##              coef exp(coef)  se(coef)      z Pr(>|z|)
## preydiv  0.957799  2.605953  6.772711  0.141    0.888
## log_sl  -0.008404  0.991631  0.083328 -0.101    0.920
## cos_ta   0.031707  1.032215  0.130643  0.243    0.808
## 
##         exp(coef) exp(-coef) lower .95 upper .95
## preydiv    2.6060     0.3837 4.477e-06 1.517e+06
## log_sl     0.9916     1.0084 8.422e-01 1.168e+00
## cos_ta     1.0322     0.9688 7.990e-01 1.333e+00
## 
## Concordance= 0.494  (se = 0.029 )
## Likelihood ratio test= 0.09  on 3 df,   p=1
## Wald test            = 0.09  on 3 df,   p=1
## Score (logrank) test = 0.09  on 3 df,   p=1

summary(ssf2)

## Call:
## coxph(formula = Surv(rep(1, 1518L), case_) ~ preydiv * log_sl + 
##     cos_ta + strata(step_id_), data = data, model = TRUE, method = "exact")
## 
##   n= 1516, number of events= 138 
##    (2 observations deleted due to missingness)
## 
##                      coef  exp(coef)   se(coef)      z Pr(>|z|)
## preydiv         2.791e+01  1.318e+12  2.177e+01  1.282    0.200
## log_sl          1.663e+00  5.273e+00  1.286e+00  1.293    0.196
## cos_ta          2.338e-02  1.024e+00  1.311e-01  0.178    0.858
## preydiv:log_sl -2.744e+00  6.433e-02  2.101e+00 -1.306    0.192
## 
##                exp(coef) exp(-coef) lower .95 upper .95
## preydiv        1.318e+12  7.589e-13 3.919e-07 4.430e+30
## log_sl         5.273e+00  1.897e-01 4.241e-01 6.556e+01
## cos_ta         1.024e+00  9.769e-01 7.917e-01 1.324e+00
## preydiv:log_sl 6.433e-02  1.555e+01 1.046e-03 3.955e+00
## 
## Concordance= 0.555  (se = 0.029 )
## Likelihood ratio test= 1.81  on 4 df,   p=0.8
## Wald test            = 1.81  on 4 df,   p=0.8
## Score (logrank) test = 1.81  on 4 df,   p=0.8

3.3 HMM

Finally, we will fit the HMM using momentuHMM (McClintock and Michelot, 2018). In preparation, we define initial parameters, and then update the parameters using our fit model, to ultimately fit a more refined model.

# prepare dataset
data_hmm <- seal %>%
  make_track(lon, lat, date) %>%
  extract_covariates(fish_raster) %>%
  mutate(ID = 1,
         x = x_,
         y = y_,
         date = t_) %>%
  dplyr::select(ID, x, y, date, preydiv) %>%
  as.data.frame()

# convert into momentuHMM format
data_hmm <- momentuHMM::prepData(data_hmm, coordNames = c("x", "y"), 
                                 type = "UTM",
                                 covNames = c("preydiv"))

# define initial parameters
nbStates <- 3 # number of states
stepDist <- "gamma" # step distribution
angleDist <- "vm" # turning angle distribution

mu0 <- c(5000, 18000, 38000) # mean step length for each state
sigma0 <- c(5000,  10000, 10000) # sd step length for each state
stepPar0 <- c(mu0, sigma0) 
kappa0 <- c(0.35, 0.55, 0.5) # turning angle for each state

formula = ~ preydiv # identify covariates

# fit basic 3-state hmm
set.seed(2023)
hmm1 <- momentuHMM::fitHMM(data=data_hmm, nbStates=nbStates,
            dist=list(step=stepDist,angle=angleDist),
            Par0=list(step=stepPar0,angle=kappa0))

# retrieve parameters to refine model
Par0_hmm1 <- momentuHMM::getPar0(hmm1, formula=formula)

# fit a refined hmm with parameters from hhm1
set.seed(2023)
hmm2 <- momentuHMM::fitHMM(data=data_hmm, nbStates=3,
            dist=list(step=stepDist,angle=angleDist),
            Par0=Par0_hmm1$Par,
            delta0 = Par0_hmm1$delta, 
            beta0 = Par0_hmm1$beta,
            formula=formula)

We plot the decoded states (estimated behaviours).

# plot decoded states
data_hmm$state <- as.factor(momentuHMM::viterbi(hmm2))
hmmstate_plot <- ggplot() + 
  scale_fill_viridis(option = "mako", limits = c(0.5, 0.75), name = "Prey\ndiversity") +
  geom_path(data=data_hmm, aes(x=x, y=y, color = state, group =ID)) + 
  geom_point(data=data_hmm, aes(x=x, y=y, color = state, shape = state), size=2, alpha = 0.8) + 
  scale_color_manual(values = c("#99DDB6", "#539D9C", "#312C66"), 
                     labels = c("Localized\nsearch", "Area-restricted\nsearch", "Travelling"), 
                     name = "HMM state") +
  scale_shape_manual(values = c(15, 16, 17), 
                     labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travelling"), 
                     name = "HMM state") +
  geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
  coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F)
hmmstate_plot

Plot a histogram of step lengths for each state. In order to do this, we need to transform the movement data to UTM and refit our HMM in order to get step length in a meaningful unit (i.e., kilometers)

## prep data in lat/lon and refit model (so step length is in kilometers)
data_hmm <- seal %>%
  make_track(lon, lat, date) %>%
  extract_covariates(fish_raster) %>%
    mutate(ID = 1,
         x = x_,
         y = y_,
         date = t_) %>%
  dplyr::select(ID, x, y, date, preydiv) %>%
  as.data.frame() %>%
  st_as_sf(coords = c("x", "y")) %>%
  st_set_crs("+proj=laea +lat_0=60 +lon_0=-85 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0") %>%
  st_transform("+proj=longlat +datum=WGS84") %>%
  mutate(long = unlist(map(geometry,1)),
          lat = unlist(map(geometry,2))) %>%
  st_drop_geometry()

# do momentuhmm which calculates step length in KM
data_hmm <- momentuHMM::prepData(data_hmm, coordNames = c("long", "lat"), type = "LL", covNames = c("preydiv"))


# define parameters
nbStates <- 3
stepDist <- "gamma" # step distribution
angleDist <- "vm" # turning angle distribution

mu0 <- c(5, 12, 38)
sigma0 <- c(3, 5, 8)
stepPar0 <- c(mu0, sigma0)
kappa0 <- c(0.35, 0.55, 0.5)

# fit HMM (with step in km)
set.seed(2023)
hmm1_km <- momentuHMM::fitHMM(data=data_hmm, nbStates=nbStates,
            dist=list(step=stepDist,angle=angleDist),
            Par0=list(step=stepPar0,angle=kappa0))

# get parameters 
Par0_hmm1_km <- momentuHMM::getPar0(hmm1_km, formula=formula)

# fit HMM with parameters from hmm1_km
set.seed(2023)
hmm2_km <- momentuHMM::fitHMM(data=data_hmm, nbStates=3,
            dist=list(step=stepDist,angle=angleDist),
            Par0=Par0_hmm1_km$Par,
            delta0 = Par0_hmm1_km$delta, 
            beta0 = Par0_hmm1_km$beta,
            formula=formula)

# add the state estimate from the HMM
data_hmm$state <- as.factor(momentuHMM::viterbi(hmm2_km))

# calculate frequencies of states
v <- momentuHMM::viterbi(hmm2_km)
stateFreq <- table(v) / length(v)

# plot colours
colours.states <- c("#99DDB6", "#539D9C", "#312C66")

# generate sequence for x axis of density functions
x <- seq(0, 50, length=1000)

# get converged mean and sd for each state 
meanARS <- hmm2_km$mle$step[1,1]  
sdARS <- hmm2_km$mle$step[2,1]    

meanCR <- hmm2_km$mle$step[1,2]    
sdCR <- hmm2_km$mle$step[2,2]     

meanTR <- hmm2_km$mle$step[1,3]    
sdTR <- hmm2_km$mle$step[2,3]    

# calculate shape and scale of the gamma distributions from mean and sd
sh <- function(mean, sd) { return(mean^2 / sd^2)}
sc <- function(mean, sd) { return(sd^2 / mean)}

# get density functions of the distributions
y_ARS <- dgamma(x, shape=sh(meanARS,sdARS),  scale=sc(meanARS,sdARS)) * stateFreq[[1]]
y_CR <- dgamma(x, shape=sh(meanCR,sdCR),  scale=sc(meanCR,sdCR)) * stateFreq[[2]]
y_TR <- dgamma(x, shape=sh(meanTR,sdTR),  scale=sc(meanTR,sdTR)) * stateFreq[[3]]


# combine densities in a single dataframe for more convenient plotting
df.y_ARS <- data.frame(dens=y_ARS, State="Foraging", x=x)
df.y_CR <- data.frame(dens=y_CR,  State="ARS", x=x)
df.y_TR <- data.frame(dens=y_TR,  State="Travelling", x=x)
statedis <- rbind(df.y_ARS, df.y_CR, df.y_TR)

# plot distributions
hmmstepdist_plot <- ggplot() +
  geom_line(data=statedis,aes(x=x,y=dens,colour=State,linetype=State), size=1.2) +
  scale_colour_manual(values=c(colours.states,"#000000"), 
              breaks = c('Foraging', 'ARS', 'Travelling'), 
              labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travelling")) +
  scale_linetype_manual(values=c("solid","solid", "solid"), 
              breaks = c('Foraging', 'ARS', 'Travelling'), 
              labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travelling")) +
  ylab("Density") + 
  xlab("Step length (kms)") +
  theme_minimal()

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.

## Warning: Please use `linewidth` instead.

hmmstepdist_plot

Plot a histogram of turning angle for each state.

# generate sequence for x axis of density functions
x <- seq(-pi, pi,length=1000)

# get converged mean and concentration for each state 
meanARS <- hmm2_km$mle$angle[1,1]  
sdARS <- hmm2_km$mle$angle[2,1]    

meanCR <- hmm2_km$mle$angle[1,2]    
sdCR <- hmm2_km$mle$angle[2,2]  

meanTR <- hmm2_km$mle$angle[1,3]    
sdTR <- hmm2_km$mle$angle[2,3]  

# get density functions of the distributions
y_ARS <- CircStats::dvm(x, mu=meanARS,  kappa=sdARS) * stateFreq[[1]]
y_CR <- CircStats::dvm(x, mu=meanCR,  kappa=sdCR) * stateFreq[[2]]
y_TR <- CircStats::dvm(x, mu=meanTR,  kappa=sdTR) * stateFreq[[3]]


# combine densities in a single dataframe for more convenient plotting
df.y_ARS <- data.frame(dens=y_ARS, State="Foraging",x=x)
df.y_CR <- data.frame(dens=y_CR,  State="ARS",x=x)
df.y_TR <- data.frame(dens=y_TR,  State="Travelling",x=x)

cmb <- rbind(df.y_TR, df.y_CR, df.y_ARS)

# plot distributions
hmmangledist_plot <- ggplot() +
  geom_line(data=cmb,aes(x=x,y=dens,colour=State), size = 1.2) +
  scale_colour_manual(values=c(colours.states), breaks = c('Foraging', 'ARS', 'Travelling'), labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travelling")) +
  scale_x_continuous(limits=c(-pi,pi))+
  ylab("Density") + 
  xlab("Turn angle (radians) ") +
  theme_minimal()
hmmangledist_plot

4 Prediction maps

Now we will estimate the utilization distribution from each model to demonstrate the inference from each model. The utilization distribution is defined as the two-dimensional relative frequency distribution of space use of an animal (Van Winkle 1975). This is a simple calculation for the RSF, where we multiply the model coefficent with the resource (prey diversity), exponentiate (since it is a logistic regression), and normalize the estimate. The calculations are more complex for the SSFs since they are conditional models that integrates the movement process. Thus, for the SSFs we calculate the steady-state utilization distribution (SSUD), which is the long-term expectation of the space-use distribution across the landscape (Signer et al. 2017). Luckily, amt has functions to estimate the SSUD.

First, we prepare a dataframe to predict on.

# prep the fish data
newfish <- fish_raster %>%
  terra::as.data.frame(xy = TRUE) %>%
  filter(x > 100000 & x < 600000 & y > -550000 & y < 0)

4.1 RSF

Now we will predict the estimated probability of use from the RSF by hand.

# grab model coefficients
modcoef <- summary(rsf1)$coef

# prediction for each cell
x <- exp(modcoef[1] + (modcoef[2] * newfish$preydiv))

# scale the data
x <- x / (1 + x)

# range fn
range01 <- function(x){(x-min(x))/(max(x)-min(x))}

# set the range from zero to one
newfish$rsf_prediction <-range01(x)

# plot
map_RSF <- ggplot() + 
  geom_tile(data = newfish, aes(x = x,y = y, fill = rsf_prediction)) +
  scale_fill_viridis(option = "mako", name = "RSF prediction", limits = c(0, 0.3)) +
  geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
  coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F) 

## Regions defined for each Polygons

map_RSF

4.2 SSF

We will first estimate the SSUDs from the simple SSF that does not allow prey diversity to affect the movement kernel. We need to refit our SSF model to include terms for the home range in order to make sure our simulated track stays within our study area.

# generate availability sample
set.seed(2023)
data_ssf <- seal %>%
  mutate(date = as.POSIXct(date)) %>%
  make_track(lon, lat, date) %>%
  steps() %>%
  random_steps() %>% 
  arrange(case_) %>%
  amt::extract_covariates(fish_raster, where = "both") %>%
  na.omit() 

# fit SSF1 model
m1 <- data_ssf |> 
  fit_clogit(case_ ~ preydiv_end +  cos(ta_) + log(sl_) + 
              # add terms for home range
               x2_ + y2_ + I(x2_^2 + y2_^2) +
               strata(step_id_))

We will now simulate a track and visually observe it.

# set starting position for the simulation
set.seed(2023)
start <- make_start((seal %>%
                       mutate(date = as.POSIXct(date)) %>%
                       make_track(lon, lat, date))
                          # random sample along the track as the start
                          [sample(1:nrow(seal), 1), ], 
                    dt_ = hours(2)) 

# generate redistribution kernel
k1 <- redistribution_kernel(m1, map = fish_raster, start = start,
                            stochastic = TRUE,
                            tolerance.outside = 0.1,
                            n.control = 1e3)

# Now simulate a path of length 1e3
n_steps = 1e3
n_steps1 = n_steps + 1
set.seed(2023)
p1 <- amt::simulate_path(x = k1, n = n_steps, start = start, verbose = TRUE)

# plot simulated track
ssf_track_1 <- fishmap +
  geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
  geom_point(data = p1, aes(x = x_,y = y_), alpha = 0.61) +
  geom_path(data = p1, aes(x = x_,y = y_))

## Regions defined for each Polygons

ssf_track_1

We can see that it mostly stays within the study area. We will use this track to estimate the SSUD, and visualize the results.

# estimate SSUD
uds_ssf1 <- tibble(rep = 1:n_steps1, 
    x_ = p1$x_, y_ = p1$y_,
    t_ = p1$t_, dt = p1$dt) |> 
  filter(!is.na(x_)) |> 
  make_track(x_, y_) |> 
  hr_kde(trast = fish_raster, which_min = "global") %>%
  hr_ud() %>% 
  terra::as.data.frame(xy = TRUE)

# plot SSUD
map_SSF1 <- ggplot() + 
  geom_tile(data = uds_ssf1, aes(x = x,y = y, fill = preydiv)) +
  scale_fill_viridis(option = "mako", name = "SSF1 prediction") +
  geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
  coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F) 

## Regions defined for each Polygons

map_SSF1

We will follow the same steps to generate a SSUD for the SSF that allows prey diversity to affect the movement kernel.

# fit SSF2 model
m2 <- data_ssf |> 
  fit_clogit(case_ ~ preydiv_end +  cos(ta_)+
               preydiv_end:log(sl_) + 
               # add terms for home range
               x2_ + y2_ + I(x2_^2 + y2_^2) +
               strata(step_id_))


# set starting position for the simulation
set.seed(2023)
start <- make_start((seal %>%
                       mutate(date = as.POSIXct(date)) %>%
                       make_track(lon, lat, date))
                          # first location of the track as the start
                          [sample(1:nrow(seal), 1), ], 
                    dt_ = hours(2)) 

# generate redistribution kernel
k2 <- redistribution_kernel(m2, map = fish_raster, start = start,
                            stochastic = TRUE, 
                            tolerance.outside = 0.3,
                            n.control = 1e3)

# Now simulate a path of length 1e3
n_steps = 1e3
n_steps1 = n_steps + 1
set.seed(2023)
p2 <- amt::simulate_path(x = k2, n = n_steps, start = start)

# plot simulated track
ssf_track_2 <- fishmap +
  geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
  geom_point(data = p2, aes(x = x_,y = y_), alpha = 0.61) +
  geom_path(data = p2, aes(x = x_,y = y_)) 

## Regions defined for each Polygons

ssf_track_2

# estimate SSUD
uds_ssf2 <- tibble(rep = 1:n_steps1, 
    x_ = p1$x_, y_ = p1$y_,
    t_ = p1$t_, dt = p1$dt) |> 
  filter(!is.na(x_)) |> 
  make_track(x_, y_) |> 
  hr_kde(trast = fish_raster, which_min = "local") %>%
  hr_ud() %>% 
  terra::as.data.frame(xy = TRUE)

# plot SSUD
map_SSF2 <- ggplot() + 
  geom_tile(data = uds_ssf2, aes(x = x,y = y, fill = preydiv)) +
  scale_fill_viridis(option = "mako", name = "SSF2 prediction") +
  geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
  coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F) 

## Regions defined for each Polygons

map_SSF2

4.3 HMM

We will first estimate the stationary state probabilities of each state based on prey diversity. This is easily done using momentuHMM’s function stationary().

# grab estimated stationary state probabilities from our fitted HMM
x <- as.data.frame(momentuHMM::stationary(hmm2, data.frame(preydiv = newfish$preydiv)))
newfish$hmm_state1 <- x$state.1
newfish$hmm_state2 <- x$state.2
newfish$hmm_state3 <- x$state.3 

# prepare data
newfish_long <- newfish %>%
  tidyr::pivot_longer(cols = hmm_state1:hmm_state3, 
               names_to = "model", values_to = "prediction") %>%
  mutate(dplyr::across(model, factor, levels=
               c("hmm_state1", "hmm_state2", "hmm_state3")))

# plot
map_HMM <- ggplot() + 
  geom_tile(data = newfish_long, aes(x = x, y = y, fill = prediction)) +
  scale_fill_viridis(option = "mako", limits = c(0,1)) +
  labs(fill = 'HMM predicted\nprobability') +
  geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), colour = "white", fill = "grey90", size = 0.3) +
  coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F) +
  facet_wrap(~model, labeller = as_labeller(c('hmm_state1' = "Localized search",
                                              'hmm_state2' = "Area-restricted search",
                                              'hmm_state3' = "Travel")))

## Regions defined for each Polygons

map_HMM

5 Line plots of prey diversity vs probability of selection

5.1 RSF

base <- newfish %>% 
  #  (note, it doesn't matter what values you pick for sl_ and ta_)
  mutate(log_sl = log(45),
         cos_ta = cos(1))

x1 <- base 

x2 <- base  %>% 
  mutate(preydiv = mean(base$preydiv))

# So how do we iterate this? We could use a for loop, we could use lapply, etc.
# I will use lapply here, with a custom function inside.

log_rss_list <- lapply(1:nrow(x1), function(i) {
  # Calculate log-RSS for that row
  xx <- log_rss(rsf1, x1[i,], x2[i,], ci = "se")
  # Return the element $df
  return(xx$df)
})

# We now have a list with an element for each row. We can combine those rows
# into a single data.frame with 'dplyr::bind_rows()'
res <- dplyr::bind_rows(log_rss_list)

# PLOT
rsf_line <- ggplot(res, aes(x = preydiv_x1, y = (log_rss))) +
  geom_line(size = 1, color = "#F8D59F",linetype = 2) +
  geom_ribbon(aes(ymin=lwr, ymax=upr, x=preydiv_x1), alpha = 0.4, fill = "#F8D59F") +
  xlab("Prey diversity") +
  ylab("log-RSS") +
  theme_minimal()

5.2 SSF

## Prediction for SSF1

x1 <- base 

x2 <- base  %>% 
  mutate(preydiv = mean(base$preydiv))

# So how do we iterate this? We could use a for loop, we could use lapply, etc.
# I will use lapply here, with a custom function inside.

log_rss_list <- lapply(1:nrow(x1), function(i) {
  # Calculate log-RSS for that row
  xx <- log_rss(ssf1, x1[i,], x2[i,], ci = "se")
  # Return the element $df
  return(xx$df)
})

# We now have a list with an element for each row. We can combine those rows
# into a single data.frame with 'dplyr::bind_rows()'
res <- dplyr::bind_rows(log_rss_list)


# PLOT
ssf_line_1 <- ggplot() +
  geom_line(data = res, aes(x = preydiv_x1, y = log_rss), size = 1, linetype = 3) +
  geom_ribbon(data = res, aes(ymin=lwr, ymax=upr, x=preydiv_x1), alpha = 0.3) +
  xlab("Prey diversity") +
  ylab("log-RSS") +
  theme_minimal()


## Prediction for SSF2
nums = seal %>%
  make_track(lon, lat, date) %>%
  steps %>%
  summarize(quants = quantile(sl_, c(0.25, 0.5, 0.75))) %>%
  pull()

x1 <- base 

results <- lapply(nums, function(i) {
x1$log_sl <- i

x2 <- base  %>% 
  mutate(preydiv = mean(base$preydiv))

# So how do we iterate this? We could use a for loop, we could use lapply, etc.
# I will use lapply here, with a custom function inside.

log_rss_list <- lapply(1:nrow(x1), function(i) {
  # Calculate log-RSS for that row
  xx <- log_rss(ssf2, x1[i,], x2[i,], ci = "se")
  # Return the element $df
  return(xx$df)
})

# We now have a list with an element for each row. We can combine those rows
# into a single data.frame with 'dplyr::bind_rows()'
res <- dplyr::bind_rows(log_rss_list)

} )
results <- dplyr::bind_rows(results) %>%
  mutate(log_sl_x1 = as.factor(round(log_sl_x1, 1)))

# PLOT
ggplot(results, aes(x = preydiv_x1, y = (log_rss))) +
  geom_line(size = 1, aes(color = log_sl_x1, group = log_sl_x1)) +
  geom_ribbon(aes(ymin=lwr, ymax=upr, x=preydiv_x1, fill = log_sl_x1, group = log_sl_x1), alpha = 0.3) +
  xlab("Prey diversity") +
  ylab("log-RSS") +
  theme_minimal()

# Plot
ssf_line_1 +
  geom_line(data = results, size = 1, aes(x = preydiv_x1, y = (log_rss), 
                                          color = log_sl_x1, group = log_sl_x1)) +
  geom_ribbon(data = results, aes(ymin=lwr, ymax=upr, x=preydiv_x1, 
                                  fill = log_sl_x1, group = log_sl_x1), alpha = 0.3)

5.3 HMM

# grab stationary probabilities
ps <- momentuHMM::plotStationary(hmm2, plotCI= TRUE, return = TRUE)

state1 <- ps$preydiv$'state 1' %>% mutate(state = 1)
state2 <- ps$preydiv$'state 2' %>% mutate(state = 2)
state3 <- ps$preydiv$'state 3' %>% mutate(state = 3)

# bind to one data frame
pdat <- rbind(state1, state2, state3) %>%
  mutate(state = as.character(state)) 

# plot
hmmprobs_plot <- ggplot() + 
  geom_line(data = pdat, aes(x = cov, y = est, color = state)) +
  geom_ribbon(data = pdat, aes(x=cov, y=est, ymax=est+se, ymin=est-se, fill = state), 
              alpha = 0.4, show.legend = TRUE) +
  ylab("Stationary state probabilties") +
  xlab("Prey diversity") +
  scale_color_manual(values = c("#99DDB6", "#539D9C", "#312C66"), name = "HMM state",
                     labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travel")) +
    scale_fill_manual(values = c("#99DDB6", "#539D9C", "#312C66"), name = "HMM state",
                     labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travel")) +
  theme_minimal()

hmmprobs_plot

6 References

Florko, K.R.N., Tai, T.C., Cheung, W.W.L., Sumaila, U.R., Ferguson, S.H., Yurkowski, D.J., Auger-Méthé, M. 2021. Predicting how climate change threatens the prey base of Arctic marine predators. Ecology Letters, 24: 2563-2575. https://doi.org/10.1111/ele.13866

Florko, K.R.N., Shuert, C.R., Cheung, W.W.L., Ferguson, S.H., Jonsen, I.D., Rosen, D.A.S., Sumaila, U.R., Tai, T.C., Yurkowski, D.J., Auger-Méthé, M. 2023. Linking movement and dive data to prey distribution models: new insights in foraging behavior and potential pitfalls of movement analyses. Movement Ecology, 11:17 https://doi.org/10.1186/s40462-023-00377-2

McClintock, B.T., Michelot, T. 2018. momentuHMM: R package for generalized hidden markov models of animal movement. Methods Ecol. Evolut. 9, 1518-1530. https://doi.org/10.1111/2041-210X.12995

Signer, J., Fieberg, J., & Avgar, T. (2017). Estimating utilization distributions from fitted step‐selection functions. Ecosphere, 8(4), e01771. https://doi.org/10.1002/ecs2.1771

Signer, J., J. Fieberg, and T. Avgar. 2019. Animal movement tools (amt): R package for managing tracking data and conducting habitat selection analyses. Ecology and Evolution 9:880–890. https://doi.org/10.1002/ece3.4823

Van Winkle, W. (1975). Comparison of several probabilistic home-range models. The Journal of wildlife management, 118-123. https://doi.org/10.2307/3800474